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PREFACE 

In  our  efforts  to  keep  abreast  of  developments  in  computa- 
tional mathematics,  we  have  found  it  useful  to  make  this  listing  of 
abstracts  of  the  many  departmental  reports  which  we  have  accumulated 
during  the  past  year. 

The  abstracts  are  largely  those  provided  with  the  documents; 
where  no  abstract  was  given  we  used  selected  sentences  from  the  intro- 
duction.  Occasionally  we  condensed  an  abstract  which  seemed  too  long 
for  our  purposes.   Sources  of  the  reports  are  identified  by  letter 
codes  prefacing  the  report  number.   A  key  to  these  codes  is  given  at 
the  end  of  the  listing. 

We  wish  to  thank  Nancy  Freece  not  only  for  her  excellent 
and  rapid  job  of  typing  but  also  for  doing  a  considerable  amount  of 
the  work  of  hunting  through  the  documents  to  assemble  the  information 
listed. 


ERROR  ANALYSIS 


Algorithms  to  Reveal  Properties   of  Floating-Point  Arithmetic 
by  Michael  A.    Malcolm 

Two   algorithms  are  presented  in  the   form  of  FORTRAN  sub- 
routines.     Each  subroutine   computes  the   radix  and  number  of  digits 
of  the   floating-point   numbers   and  whether  rounding  or  chopping  is 
done  by  the  machine  on  which   it   is   run.      The  methods   are   shown  to 
work  on  any   "reasonable"   floating-point    computer. 

STAN-CS-T1-211  March  1971  8  Pages  Program  Incl. 


On  the  Precision  Attainable  with  Various   Floating-Point   Number  Systems 
by  R.    P.    Brent 

For  scientific   computations  on   a  digital   computer  the  set 
of  real   numbers    is   usually  approximated  by  a  finite  set  F  of  "floating- 
point  numbers".      We   compare  the  numerical   accuracy  possible  with 
different   choices  of  F  having  approximately  the   same   range  and  requir- 
ing the   same  wordlength.      In  particular,   we   compare  different   choices 
of  base    (or  radix)  with  the  usual   floating-point   systems.      The  emphasis 
is  on  the   choice  of  F,    not   on  the   details  of  the  number  representation 
or  the  arithmetic,   but  both  rounded  and  truncated  arithmetic  are   consi- 
dered.     Theoretical  results   are  given,    and  some   simulations  of  typical 
floating-point   computations    (forming  sums,    solving  systems   of  linear 
equations,    finding  eigenvalues)    are   described.      If  the  leading  fraction 
bit   of  a  normalized  base-2  number  is  not   stored  explicitly   (saving  a 
bit),    and  the   criterion   is  to  minimize  the  mean  square  roundoff  error, 
then  base  2   is  best.      If  unnormalized  numbers   are  allowed,    so  the   first 
bit  must  be   stored  explicitly,    then  base  h    (or  sometimes  base  8)    is  the 
best   of  the  usual   systems. 

YORK-RC   3751  February  1972  28  Pages 


Performance  Statistics    for  the  FORTRAN  IV    (H)    and  PL/1    (Version  5) 
Libraries   in   IBM  OS/360   Release  18 
by   Kenneth   E.    Hills trom 

The   computational   subroutine  libraries   associated  with  the 
FORTRAN  IV   (H)    and  PL/I    (F)    (Version   5)    compilers   of  OS  Release  18 
have  been  tested.      The  testing  techniques   are   described  and  accuracy 
and  timing  statistics   are  presented;   proper  operation  under  error 
conditions  was   verified. 

ANL-7666  August  1970  70  Pages  Program  Incl. 


FUNCTION  EVALUATION  AND  COMPLEXITY  THEORY 


On  The  Number  of  Multiplications  for  the  Evaluation  of  a  Polynomial 
and  Some  of  its  Derivatives 
by  Mary  Shaw  and  J.  F.  Traub 

Some  recent  work  in  computational  complexity  has  dealt  with 
the  number  of  arithmetic  operations  needed  to  evaluate  a  polynomial 
or  a  polynomial  and  its  first  derivative.   Here  we  consider  the 
evaluation  of  a  polynomial  and  its  first  m  derivatives  and,  in  parti- 
cular, the  calculation  of  all  the  derivatives. 

CMU  -  No  #  August  1972  18  Pages 


Computational  Complexity  of  Iterative  Processes 
by  J.  F.  Traub 

The  theory  of  optimal  algorithmic  processes  is  part  of 
computational  complexity.   This  paper  deals  with  analytic  computa- 
tional complexity.   The  relation  between  the  goodness  of  an  iteration 
algorithm  and  its  new  function  evaluation  and  memory  requirements 
are  analyzed.  A  new  conjecture  is  stated. 

CMU-CS-71-105        October  1971  25  Pages 


On  the  Additions  Necessary  to  Compute  Certain  Functions 
by  David  G.  Kirkpatrick 

We  describe  low  level  or  functional  complexity,  fitting  it 
into  the  general  framework  of  computational  complexity  and  presenting 
a  number  of  results  in  arithmetic  complexity.  We  put  forward  a  new 
notion  of  independence,  called  rational  independence,  which  applies 
to  the  terms  of  a  computed  expression.  We  are  able  to  show  that  this 
notion  serves  as  a  measure  for  the  number  of  additions  needed  to 
compute  a  function. 

TOR-Tech.  Rpt .  39         February  1972       75  Pages     [M.S.  Thesis] 


Some  Results  in  the  Study  of  Algorithms 
by  J.  Ian  Munro 

Functions  such  as  polynomials,  matrix  products,  and  digraph 
transitive  closures  are  in  practice  computed  very  often.   In  recent 
years,  considerable  attention  has  been  turned  to  efficient  computation 
of  such  functions,  with  the  hope  of  minimizing  the  time  needed  to 
compute  them  on  a  general  purpose  digital  computer  as  we  now  know  it. 


In  this  thesis  several  new  algorithms  are  given  for  the  evaluation  of 
such  functions,  and  in  some  cases  these  are  shown  to  be  optimal  under 
a  reasonable  model  of  computation. 

It    is   shown  that   if  a  polynomial,   or  set   of  polynomials,    is 
to  be  evaluated  at   a  large  number  of  points,    and  these  points   are  all 
given  at  the   same  time,   then  the  number  of  arithmetic  operations 
needed  per  "polynomial-point"   is   less  than  linear  in  the   degree  of 
the  polynomial. 

An  algorithm  requiring  work  of  the  order  of  max   (n  log  n,e) 
is   given   for  the   determination  of  the   strongly  connected  components  of 
an  n-node,    e-edge   digraph. 

Attention   is  turned  to  problems  of  computing  arithmetic 
functions   under  a  model  of  -  computation   in  which  a  large  degree  of 
parallelism  is  permitted.      New  optimal   and  almost   optimal   algorithms 
are  presented  for   several   such   computations. 

TOR-Tech.    Rpt .    32  October  1971  102  Pages  Ihesis 


Application  of  Continued  Fractions    for  Fast   Evaluation  of  Certain 
Functions   on   a  Digital   Computer 
by  Amnon   Bracha 

The  purpose  of  this   paper  is   to    develop  a  method  for  evalua- 
tion of  certain  elementary  functions  on  a  digital   computer  by  the  use 
of  continued  fractions.      The  time  required  for  this   evaluation  is   drasti- 
cally reduced  by  using   "short"  operations   like   shift  and  add,   instead  of 
multiplications.      Functional   consistency  is  the  most    important   factor 
that    allows   the  expansion  of  a  function   into   a  continued  fraction. 
Several   cases   are   discussed  and  in  particular  the  solution  of  the   quadra- 
tic  equation   is    discussed  in  more   detail  to   demonstrate  the   convergence 
of  the  method. 

UIUC-R- 72-510  March  1972  29  Pages 


Bounds   on  Polynomial   Evaluation  Algorithms 
by  Larry  Joseph  Stockmeyer 

The  purpose  of  this  work  is  to   investigate  the  number  of 
arithmetic   operations   required  by  algorithms  which  evaluate   polynomials. 
Previous   results   show  that    a  polynomial  of  degree  n   requires   at   least 
n/2  multiplication/divisions   and  at   least   n  addition/subtractions    for 
its   evaluation   if  the  coefficients   of  the  polynomial   are   suitably 
independent    irrational  numbers.      However,   the  coefficients   of  any  poly- 
nomial that  would  be   evaluated  in  practice  are  represented  only  to  a 
finite  accuracy  and  are  therefore  rational  numbers.      The   above  results 


are   extended  to   show  that   the  same  lower  bounds  hold  for  almost   all 
rational  polynomials    if  the  polynomial   is  being  evaluated  efficiently. 
Another  lower  bound  result    is   given  that   shows  that   almost   all   rational 
polynomials   of  degree  n  require  at   least    vn  multiplication/divisions 
for  their  evaluation  by  any  algorithm,    efficient   or  not. 

Several   algorithms   are  presented  which   can   in  theory 
evaluate  any  rational  polynomial  using  0(/n)   multiplications   and  many 
additions.      While  of  no  practical  use   for  rational   polynomials    in 
general,   these  algorithms   do  turn  out  to   give  methods   for  evaluating 
a  polynomial   at   a  matrix  argument  which   are  more  efficient  than  pre- 
vious methods. 

MAC   TR-98  April  1972  5*+  Pages 


APPROXIMATION 


Variational  Study  of  Nonlinear  Spline   Curves 
by  E.    H.    Lee  and  G.    E.    Forsythe 

This   is   an  exposition  of  the  variational   and  differential 
properties   of  nonlinear  spline   curves,  based  on  the  Euler-Bernoulli 
theory  for  the  bending  of  thin  beams  or  elastica.      For  both  open   and 
closed  splines  through  prescribed  nodal  points   in  the  euclidean  plane, 
various   types   of  nodal   constraints    are  considered,    and  the   correspond- 
ing algebraic   and  differential   equations   relating  curvature,   angle, 
arc  length,    and  tangential   force   are   derived  in   a  simple  manner. 
The  results    for  closed  splines   are  apparently  new,    and  they   cannot 
be  derived  by  the  consideration  of  a  constrained  conservative  system. 
There   is   a  survey  of  the  scanty  recent   literature. 

STAN-CS-229-71  August  1971  26  Pages 


An  Algorithm  for  Fitting  Related  Sets   of  Straight-Line  Data 
by   Geneva  G.    Bel ford 

Many  physical   experiments   give  rise  to   sets   of  curves   related 
by  the   requirement  that,    although   certain  of  the  curve  parameters  may 
vary  from  curve  to   curve,    others   should  be  the   same   for  all  of  the 
curves.      To  get   the   "best"   values   of  the  common  parameters,    one  would 
like  to    fit    all   of  the   curves   simultaneously  by  the  appropriate  theore- 
tical  expressions.      This   paper  deals  with  this   problem,   presenting  an 
algorithm  for   its   solution,    in  the  case  that  the   curves   are   straight 
lines  with   common   slope  and  "best"   fit    is    defined  in  the  uniform   (or 
minimax)    sense. 

CAC   No.    5U  December  1972  12  Pages 


An  Implementation  of  the  Remes  Algorithm  for  Minimax  Approximation 
by  Roland  Olofsson 

The  ALGOL  procedure   for  the  Remes    algorithm  given  by  G.    H. 
Golub   and  L.    B.    Smith    (Comm  ACM  Ik    (1971 ),    737-7^6  Alg  klk)   has  been 
implemented  and  tested  on   a  CD  3200   computer.      Some  examples   of  ill 
conditioned  generalized  polynomials  have  been  used  to  test  the  algorithm. 
A  comparison   is   also  made  to   another  program  based  on  the  Polya  algorithm. 

UMlNF-20.72  April   1972  Ik   Pages 


On  the  Existence  and  Characterization  of  Minimal  Projections 
by  P.    D.    Morris   and  E.    W.    Cheney 

A  projection  from  a  normed  linear  space  X  onto  a  subspace  Y 
is   a  bounded  linear  operation  P:      X  -*■  Y  having  the  property  that 
Py  =  y  for  all  y  e  Y.      Projections   play  an   important   role  in  numerical 
analysis   and  approximation  theory.      The  use  of  projections   is  based 
upon  the  acceptance  of  Px  as   an  approximation  of  x  in  the  subspace  Y. 
The  quality  of  approximations  produced  by  P  depends  upon  ||p||  and    ||I-P||, 
The  present  work  is   devoted  mainly  to  the   characterization  problem  of 
minimal  projections. 

CNA-37  January  1972  32  Pages 


Stability  Properties  of  Trigonometric   Interpolating  Operators 
by  P.    D.   Morris   and  E.   W.    Cheney 

We  consider  the  familiar  process  of  trigonometric  interpola- 
tion with  2n+l  equally-spaced  nodes.      This   process   is   interpreted  here 
as   a  linear  projection  operator  P  acting  on  the   space  C  of  all  real, 
continuous,    27r-periodic   functions.      The  range  of  this  projection  P  is 
the   subspace  n  of  all  trigonometric  polynomials  of  order  n.      Many  other 
useful  projections   exist  which  map  C  onto     II,    and  attention  here   is 
focused  on  their  extremal  properties. 

CNA-51  August  1972  20  Pages 


Extremal  Properties  of  Approximation  Operators 
by  K.    H.    Price  and  E.    W.    Cheney 

In  1910,    de  La  Vallee  Poussin  published  some  researches  on 
the   following  approximation  problem:      given  the  continuous   function 
x  =  x(t)    defined  on  the   interval    [-1,1],    and  given  n+1  points 
t    ,...,t      in   [-1,1],    find  the  polynomial  p  of  degree  <  n  which 
minimizes     max    |x(t.)   -  p(t.)|. 

0<i<n    X  X 

The  polynomial  p  depends  linearly  on  x.  In  this  respect, 
the  process  resembles  Lagrange  interpolation;  both  processes  define 
linear  projection  operators. 

In  this  paper  we  study  operators  of  this  type,  which  depend 
upon  n+1  quanta  of  information  about  the  function  x.   We  examine  the 
family  of  all  such  operators,  and  seek  out  those  that  are  optimal  in 
some  sense. 

CNA-51^  August  1972  22  Pages 


A  Generalization  of  the  Construction  Method  of  Laurent   and  Carasso 
for  Interpolating  L-Splines 
by  M.    J.    Munteanu 

The  purpose  of  this   study  is  to  generalize  the  method  given 
by  P.    T.    Laurent   and  C.    Carasso   for  the  construction  of  ordinary  poly- 
nomial  spline   interpolation  functions. 

We   explain   in  Section  2   in   detail  the  extension  of  the  method 
to  the   case  of  polynomial  spline   functions   interpolating  in  the  sense 
of  Hermite.      In  Section   3  we  treat  briefly  the   generalization  to 
interpolating  L-splines.      Finally,    in  Section  h,  we   give   a  numerical 
example  to   illustrate   certain  aspects   of  the  method. 

CNA-11  February  1971  17  Pages 


Computation  of  Smoothing  and  Interpolating  Natural  Splines   via  Local 

Bases 

by  Tom  Lyche  and  Larry  L.    Schumaker 

The  purpose  of  this  paper  is  twofold.      First  we  show  how 
smoothing  splines   can  be  represented  in  terms  of  local   support  basis 
splines,    and  how  the  corresponding  coefficients   can  be   computed  by 
solving  a  linear  system  whose  matrix  is  banded.      This  leads  to   an 
efficient   and  accurate  algorithm  for  computing  smoothing  splines. 
Secondly,    we   extend  some   ideas   of  de  Boor  to  obtain  recursion  relations 
for  the  basis   functions  which   allows  their  computation  as  positive 
combinations  of  positive  quantities. 

CNA-17  April  1971  21  Pages 


On  a  Method  of  Carasso   and  Laurent    for  Constructing  Interpolating  Splines 
by  M.    J.    Munteanu  and  L.    L.    Schumaker 

Carasso  and  Laurent   studied  a  method  for  computing  natural 
polynomial   splines   interpolating  simple   data.      We   discuss   several   similar 
methods  which  can  be   applied  to  numerical   construction  of  more  general 
interpolating  splines,    including  Lg-splines   interpolating  Extended- 
Hermit  e-Birkhoff  data. 

CNA-^0  April  1972  20  Pages 


Direct   and  Inverse  Theorems    for  Multidimensional  Spline  Approximation 
by  M.    J.    Munteanu  and  L.    L.    Schumaker 

A  general  approximation  theoretic  result   of  Butzer  and  Scherer 
is   used  to  obtain   direct   and  inverse  theorems   as  well   as   saturation 
theorems    for  multidimensional   spline  approximations   of  functions    in  cer- 
tain Sobolev  spaces. 

CNA-I13  April  1972  l6  Pages 


LINEAR  ALGEBRA 


On  the  Matrix  Polynomial,    Lambda-Matrix  and  Block  Eigenvalue  Problems 
by  J.    E.    Dennis,   Jr.,   J.    F.    Traub,    and  R.    P.   Weber 

A  matrix  S   is   a  solvent   of  the  matrix  polynomial 

M(X)    =  X111  +  A^'1  +    ...    +  A   ,    if  M(S)   =  0,   where  A.,   X  and  S   are 
1  mi 

square  matrices.      We  present   some  new  mathematical  results   for  matrix 

polynomials,    as  well  as   a  globally   convergent   algorithm  for  calculating 

such  solvents. 

In  the  theoretical  part   of  this  paper,    existence  theorems   for 
solvents,    a  generalized  division,    interpolation,   a  block  Vandermonde, 
and  a  generalized  Lagrangian  basis   are   studied.      Algorithms   are  presented 
which   generalize  Traub's   scalar  polynomial  methods,    Bernoulli's  method, 
and  eigenvector  powering. 

The  related  lambda-matrix  problem,   that   of  finding  a  scalar   \ 

such  that  . 

IAm+AAm"1+    ...    +  A 
1  m 

is   singular,    is   examined  along  with  the  matrix  polynomial   problem. 

The  matrix  polynomial   problem  can  be  cast    into   a  block  eigen- 
value  formulation  as   follows.      Given   a  matrix  A  of  order  mn,   find  a 
matrix  S  of  order  n,    such  that   AV  =  VX,   where  V  is   a  matrix  of  full 
rank.      Some  of  the   implications   of  this   new  block  eigenvalue   formulation 
are   considered. 

CMU-CS-T1-110  December  1972  117  Pages  Programs   Incl . 


The  QL-Algorithm  for  computing  the   eigenvalues  of  a  skewsymmetric  matrix 
by  Axel  Ruhe 

Two  ALGOL  procedures   are  presented  which   compute  the   eigen- 
values  and,    if  desired,   the   eigenvectors   of  a  real    skewsymmetric  matrix. 
The   first   reduces   the  matrix  into  tridiagonal    form  by  the  method  of 
Householder,    and  the   second  uses   the   QL-algorithm  to   compute  the  eigen- 
values of  the  tridiagonal  matrix.      The  procedures   take  advantage  of  the 
skewsymmetry,    and  work   considerably  faster  than   a  code  applicable  to 
general  matrices. 

UMIWE-16.72  January  1972  17  Pages  Programs   Incl. 


The  QR  algorithm  to   find  the  eigensystem  of  a  skewsymmetric  matrix 
by  Bo   Kagstrom 

The  QR   algorithm  is    specialized  in  order  to   deal   effectively 
with  a  real  skewsymmetric  matrix.      The   algorithm  starts  with  Householder 


10 


tridiagonalization,   and  uses   a  double   shift   QR  method  to  get   the 
purely  imaginary  eigenvalues,    entirely  using  real   arithmetic.      Differ- 
ent strategies   are  compared,   on  a  couple  of  numerical   examples. 

UMINF-lU.71  December  1971  lh   Pages 


Multiple   step  gradient    iterative  methods    for  computing  eigenvalues 
of  large  symmetric  matrices 
by  Britt   Marie  Alsen 

Solving  the   eigenproblem  for  the  matrix  A,    can  be  thought 
of  as   extremizing  the   Rayleigh  quotient. 

■of    \   -  x  ^ 
R(x)   =  — —  . 

X    X 

Application  of  an  optimization  algorithm  to   R(x)    finds  one   of  the   end 
eigenvalues.      In  the  present   contribution  algorithms   are  tested,  which 
extremize   R(x)    simultaneously  over  m  orthogonal  vectors,   thus   finding 
m  eigenvalues.      A  version  of  the  conjugate  gradient   algorithm  and 
steepest    descent   are  tried  on   a  couple  of  numerical   examples,   mainly 
such  emanating  from  difference  approximations. 

UMINF-15.71  December  1971  19  Pages 


Simultaneous   inverse  iteration  for  computing  orthogonal   eigenvectors 
by  Anders   Isacsson 

Inverse  iteration  is  an   accurate  and  fast  method  to   find 
eigenvectors   of  a  symmetric  matrix,   where  the  eigenvalues   are   computed 
by  bisection  or  the   QR-algorithm.      If  the  matrix  has   a  cluster  of  eigen- 
values,  the  vectors   computed  are  however  not   orthogonal.      In  the   present 
contribution   an   algorithm  is   developed  which   iterates   simultaneously  on 
a  set   of  vectors,    and  thereafter  orthogonalizes  them.      This    algorithm 
produces  orthogonal   eigenvectors  to   a  real   symmetric  matrix  having  an 
arbitrary  distribution  of  eigenvalues. 

UMIKF-19.72  March  1972  20   Pages 


Quasi -Newt on   algorithms   applied  to  eigenvalue  problems 
by  Johnny  Widen 

In  the  present   contribution  the  eigenproblem   (linear  or 
nonlinear),    (l)   T(A)x  =  0,    is   solved  by  seeking  local  minima  of  the 
functional,    (2)    HtCAJxIL  by  means   of  quasi-Newton  algorithms. 

The   experimental  evidence  suggests  that   this    is  not   a   feasible 
approach  to  the  problem,    since  even   in  well   conditioned  cases,  the  mini- 
mization of   (2)   is   difficult  numerically,   because   several  values  of  X  and 
vectors  x  give  the   function    (2)   values  of  equally  small  magnitudes. 

UMINF-23.72  April  1972  25  Pages 
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Simultaneous   iteration   for  finding  eigenvalues  of  a  symmetric  matrix 
by  Lennart   Lindstrom 

The  algorithm  of  simultaneous   iteration   as   described  "by 
Bauer  and  Rutishauser  is  tested  numerically  on  some  test  matrices. 

UMINF-13.71  November  1971  10  Pages 


MIDAS  -   Solution  of  Linear  Algebraic   Equations 
by  P.    A.    Businger 

The   information  contained  in  this   abstract    should  be  suffi- 
cient to  use  the  MIDAS  package  to   solve  general  nonsingular  systems  of 
linear  algebraic   equations,    invert  matrices,    and  compute  determinants. 
("MIDAS"  stands    for  Matrix,    Inverse,    Determinant,    and  Algebraic   System), 
The   computational  method,    nonstandard  calling  sequences,    and  an   evalua- 
tion of  the  program  are  presented  in  the  memorandum.      Error  bounds   on 
the  solution  or  inverse   are  available  as    an  option. 

BTL  No   #  January  1970  Ul   Pages  Programs   Incl 


The  Rotation  of  Eigenvectors  by  a  Perturbation  -  III 
by  Chandler  Davis  and  W.  M.  Kahan 

When  a  hermit ian  linear  operator  is  slightly  perturbed,  by 
how  much  can  its  invariant  subspaces  change?   Given  some  approximations 
to  a  cluster  of  neighbouring  eigenvalues  and  to  the  corresponding  eigen- 
vectors of  a  real  symmetric  matrix,  and  given  an  estimate  for  the  gap 
that  separates  the  cluster  from  all  other  eigenvalues,  how  much  can  the 
subspace  spanned  by  the  eigenvectors  differ  from  the  subspace  spanned  by 
our  approximations?   These  questions  are  closely  related;  both  are 
investigated  here.   The  difference  between  the  two  subspaces  is  charac- 
terized in  terms  of  certain  angles  through  which  one  subspace  must  be 
rotated  in  order  most  directly  to  reach  the  other. 

TOR  Tech.  Rpt .  6         November  1968  103  Pages 


Inclusion  Regions  for  Partitioned  Matrices 
by  Douglas  Dale  Olesky 

A  region  in  the  complex  plane  containing  all  of  the  eigen- 
values of  a  complex  matrix  A  is  said  to  be  an  inclusion  region.  The 
object  of  this  thesis  is  to  define  algorithms  for  the  computation  of 
inclusion  regions  by  partitioning  the  matrix  A. 

TOR  Tech.  Rpt.  31        September  1971      Thesis      117  Pages 
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Square  Roots   of  an  Orthogonal  Matrix 
by  Erold  ¥.    Hinds 

This  thesis   is   concerned  with  the  problem  of  determining 
numerically  the  positive   stable  orthogonal    square  root   Q  of  a  given 
real  rotation  A,   that   is,    a  rotation  Q  such  that   Of  =  A  and  Q+Q'    is 
positive   semi-definite.      Detailed  accounts  of  two  particular  algorithms 
which  may  be  used  to   construct   such  half-rotations   are  presented  herein. 

TOR  Tech.    Rpt .    37  December  1971  Thesis  167  Pages  Programs   Incl 


On  a  Characterization  of  the  Best    I     Scaling  of  a  Matrix 
by   G.    H.    Golub   and  J.    M.    Varah 

This   paper  is    concerned  with  best   two-sided  scaling  of  a 
general   square  matrix,    and  in  particular  with  a  certain   characterization 
of  that   best    scaling:      namely  that  the   first   and  last    singular  vectors 
(on  left    and  right)   of  the   scaled  matrix  have   components   of  equal 
modulus.      Necessity,    sufficiency,    and  its   relation  with  other  charac- 
terizations  are   discussed.      Then  the  problem  of  best   scaling  for  rectan- 
gular matrices   is   introduced  and  a  conjecture  made  regarding  a  possible 
best   scaling.      The  conjecture   is   verified  for  some  special   cases. 

UBC  Tech.    Rpt.    72-08  12   Pages 


On  the   Solution  of  Block  Tridiagonal   Systems  Arising  from  Certain 
Finite-Difference  Equations 
by  J.    M.    Varah 

UBC  Tech.    Rpt.    72-02  March  1972  17  Pages 


On  the  Numerical  Solution  of  Ill-Conditioned  Linear  Systems  with 
Applications  to   Ill-Posed  Problems 
by  J.    M.    Varah 

UBC  Tech.    Rpt.    71-02  October  1971  21  Pages 


Parallel   Computation  of  Eigenvalues   of  Real   Matrices 
by  David  J.    Kuck  and  Ahmed  Sameh 

This   paper  describes  the   implementation  of  three  standard 
matrix  eigenvalue  computation  methods   on  an   array  machine  with  high 
efficiency.      A  brief  description  of  the   ILLIAC   IV  computer  is  provided 
as  background  material.      Three  major  sections    follow — the   first   two 
describe  Jacobi   and  Householder  algorithms    for  real   symmetric  matrices, 
and  the  third  describes  the  QR  algorithm  for  real  nonsymmetric  matrices, 
Each  of  these  sections   is   divided  into   four  parts.      The  theoretical 
background  of  the  method  is   presented  first.     An  ILLIAC  IV  implementa- 
tion of  the  algorithm  is  presented  and  timing  estimates   are  included. 
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Next,   the   efficiency   (ratio  of  number  of  computations   actually- 
performed  to  number  of  computations   possible   in  a  given  time)   of 
the  computation  on  an  array  machine   is   derived  for  primary  memory 
constrained  matrices.      Finally,    eigenvalue   computations   for  matrices 
too  large  to  be   contained  in  primary  memory  are   discussed  in  terms 
of  a  head-per-track  secondary  storage   device. 

It    is   shown  that   for  Jacobi   and  Householder  algorithms, 
parallel   computation   efficiencies   of  90%  are   possible,   while  those 
of  the  QR  algorithm  are  rather  low. 

CAC  No.    9  November  1971  39  Pages 


Calculating  the  Eigenvectors   of  Diagonally  Dominant  Matrices 
by  M.    M.    Blevins   and  G.    W.    Stewart 

This   paper  proposes   an  algorithm  for  calculating  the  eigen- 
vectors of  a  diagonally   dominant  matrix  all  of  whose   elements    are  known 
to  high  relative  accuracy.      Eigenvectors   corresponding  to  pathologi- 
cally  close  eigenvalues   are  treated  by  computing  the  invariant   subspace 
that  they   span.      If  the  off-diagonal   elements   of  the  matrix  are   suffi- 
ciently small,   the  method  is    superior  to   standard  techniques,    and  indeed 
it  may  produce  a  complete  set   of  eigenvectors  with   an  amount   of  work 
proportional   to  the   square  of  the  order  of  the  matrix.      An  analysis   of 
the  effects  of  perturbations   in  the  matrix  on  the  eigenvectors    is   given. 

CNA-i+7  June  1972  25   Pages 


Error  and  Perturbation  Bounds   for  Subspaces   Associated  with  Certain 
Eigenvalue  Problems 
by   G.    W.    Stewart 

This   paper  describes   a  technique   for  obtaining  error  bounds 
for  certain  characteristic   subspaces   associated  with  the  algebraic 
eigenvalue  problem,    the   generalized  eigenvalue  problem,    and  the   singular 
value   decomposition.      The  method  also   gives   perturbation  bounds    for 
isolated  eigenvalues   and  useful   information   about    clusters   of  eigen- 
values.     The  bounds   are  obtained  from  an   iterative  process   for  generat- 
ing the   subspaces    in   question,    and  one  or  more   steps   of  the   iteration 
can  be  used  to  construct  perturbation  estimates  whose  error  can  be  bounded, 

CNA-1+6  June  1972  3h  Pages 


A  Note  on  Non-Hermit ian  Perturbations   of  Hermit ian  Matrices 
by  G.    W.    Stewart 

Let  A  and  E  be  Hermitian  matrices.      It    is  well  known  that  the 
eigenvalues    a.    of  A   and  8.    of  A  +  E  can  be  put   into   a  one-one   corres- 

ii  2  i2 

pondence   in   such   a  way  that    |a.    -   8.  |    <    ||e|I     and  J(a.    -   6.)      <    ||e||      . 
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In  this  note  similar  results  are  established  for  a  class  of  non- 
Hermit  ian  perturbat  ions . 

CNA-Ul  April   1972  k  Pages 


Conjugate  Direction  Algorithms   and  Matrix  Decompositions 
by  G.    W.    Stewart 

The  notion  of  a  set   of  directions   conjugate  with  respect  to 
a  matrix  A  is   generalized  and  applied  to  the   solution  of  the  linear 
system  Ax  =  b.      Variants   of  an  algorithm  for  producing  conjugate   direc- 
tions  are   shown  to  yield  several   standard  matrix  decompositions,    such 
as   the  LU  and  QR  factorizations   of  A,    as  well   as    a  number  of  previously 
studied  iterations    for   solving  linear  systems.      The   application  of  these 
results   to   systems   of  nonlinear  equations   and  to  least    squares   problems 
is   discussed. 

CNA-1+9  July  1972  16  Pages 


Linear  Stationary  Second-Degree  Methods    for  the  Solution  of  Large 

Linear  Systems 

by  David  M.    Young  and  David  R.    Kincaid 

In  this   paper  we   consider  a  class   of  linear  stationary   second- 
degree  methods   for  solving  the  linear  system  Au  =  b .      The  second-degree 
methods  which  we   shall   consider  are  related  to  ,aN"basic"linear  stationary 
method  of  first   degree  of  the   form  u^n  =   Gu         +  k  where   G  is   a  matrix 

such  that   I-G  is  nonsingular  and  k  =    (l-G)A~-Lb.      These   second-degree 
methods   q^q  0f  the   form 

u^   =  PlGU<°>   ♦  V<°>   ♦  v(l> 

(n+1)  _    (n+  (n)  (n-l)    .  .      0 

u  =  pGu  +  qu  +   ru  +  w,        n  =  1,   2,    ....,    . 

CNA-52  August   1972  29  Pages 


Error  Bounds    for  Approximate  Subspaces   of  Closed  Linear  Operators   in 
Hilbert   Space 
by  G.    W.    Stewart 

It    is  well  known  to   specialists   in  matrix  computations  that 
the   eigenvectors  of  a  matrix  corresponding  to   a  set   of  poorly  separated 
eigenvalues   are  quite   sensitive  to  perturbations    in  the  elements   of  the 
matrix.      It    is    also  known  that   for  Hermitian  matrices   the   invariant    sub- 
space  corresponding  to   a  cluster  of  eigenvalues   is   insensitive  to   such 
perturbations.      It    is  the  object   of  this   paper  to   extend  such  results   to 
nonnormal  matrices   and,   more   generally,   to   closed  operators   in  Hilbert    space 

CNA-1+  October  1970  25  Pages 
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On  the  Numerical  Properties  of  an  Iteration  for  Computing  the  Moore- 
Penrose  Generalized  Inverse 
by  G.   W.    Stewart 

In  this   paper  some  of  the  numerical   problems   associated 
with  computing  the   generalized  inverse  of  a  matrix  are  discussed  and 
illustrated  by  a  detailed  analysis  of  an   iteration  of  Ben-Israel   and 
Cohen. 

CNA-12  March  1971  27  Pages 


On  the  Solution  of  Large  Systems   of  Linear  Algebraic  Equations  with 
Sparse,    Positive  Definite  Matrices 
by  David  M.    Young 

The  object   of  this   paper   is  to   discuss  the   current   status 
of  iterative  methods   for  solving  large  systems  of  linear  algebraic 
equations.      Primary  emphasis   is  on  those  systems    involving  sparse 
matrices  where  iterative  methods   appear  more   attractive  than   direct 
methods. 

CNA-55  August  1972  U5  Pages 


The   Solution  of  Large  Systems  of  Linear  Algebraic   Equations  by 
Iterative  Methods 
by  David  M.   Young 

In  this   paper  we  provide   a   summary  of  some  of  the   available 
iterative  methods    for  solving  large   systems   of  linear  algebraic   equa- 
tions.     We   are  primarily  concerned  with  the   case  where  the  matrix  of 
the   system  is   "sparse",    i.e.,   has   only  a  few  non-zero   elements   in  each 
row.      Systems   of  this  kind  typically  arise   in  the   solution  by   finite 
difference  methods   of  elliptic  partial   differential   equations: 

CNA-19  May  1971  23  Pages 


An  Algorithm  for  the   Generalized  Matrix  Eigenvalue  Problem  Ax  =   ABx 
by  C.    B.    Moler   (Stanford  University)    and  G.   W.    Stewart 

A  new  method,    called  the  QZ   algorithm,    is  presented  for  the 
solution   of  the  matrix  eigenvalue  problem  Ax  =   ABx  with  general   square 
matrices  A  and  B.      Particular  attention  is   paid  to  the   degeneracies 
which  result  when  B  is    singular.      No  inversions   of  B  or   its   submatrices 
are   used.      The  algorithm  is    a  generalization  of  the  QR   algorithm,    and 
reduces  to   it  when  B  =   I.      A  Fortran  program  and   some   illustrative 
examples   are  included. 

CNA-32  October  1971  50  Pages 

and  STAN-CS-2U2-71  Programs   Incl. 
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A  Generalization  of  the  LR  Algorithm  to  Solve  Ax  =   ABx 

by  Linda  Kaufman 

In  this  paper,   we  will  present   and  analyze  an  algorithm  for 
finding  x  and   A   such  that  Ax  =   ABx  where  A  and  B  are  n  x  n  matrices. 
The   algorithm  does  not   require  matrix  inversion,    and  may  he  used  when 
either  or  "both  matrices   are   singular.     Our  method  is   a  generalization 
of  Rutishauser 's   LR  method  for  the   standard  eigenvalue  problem  Ax  =   Ax 
and  closely  resembles  the   QZ   algorithm  given  by  Moler  and  Stewart    for 
the   generalized  problem  given  above.      Unlike  the  QZ  algorithm,   which 
uses   orthogonal  transformations,   our  method,   the  LZ   algorithm,   uses 
elementary  transformations.      When  either  A  or  B   is    complex,   our  method 
should  be  more   efficient. 

STAN-CS- 72-276  April  1972  72  Pages  Program  Incl. 


Richardson's  Non- Stationary  Matrix  Iterative  Procedure 
by  R.    S.    Anderssen  and  G.    H.    Golub 

Because  of  its   simplicity,    Richardson's  non-stationary   itera- 
tive scheme  is   a  potentially  powerful  method  for  the  solution  of   (linear) 
operator  equations.      However,    its   general   application  has   more  or  less 
been  blocked  by    (a)      the  problem  of  constructing  polynomials,  which 
deviate  least   from  zero   on  the  spectrum  of  the  given  operator,    and  which 
are  required  for  the  determination  of  the   iteration  parameters  of  the 
non-stationary  method,    and   (b)     the   instability  of  this   scheme  with 
respect   to  rounding  error  effects. 

STAN-CS- 72-30U  August   1972  73  Pages 


Bi  diagonal  iz  at  ion  of  Matrices    and  Solution  of  Linear  Equations 
by  C.    C.    Paige 

An  algorithm  given  by  Golub   and  Kahan   for  reducing  a  general 
matrix  to  bidiagonal   form  is   shown  to  be  very  important    for  large   sparse 
matrices.      The  singular  values  of  the  matrix  are  those  of  the  bidiagonal 
form,    and  these   can  be  easily   computed.      The  bidiagonal  iz  at  ion   algorithm 
is   shown  to  be  the  basis   of  important  methods    for  solving  the  linear  least 
squares   problem  for  large  sparse  matrices.      Eigenvalues  of  certain  2-cyclic 
matrices    can   also  be   efficiently   computed  using  this  bidiagonalization. 

STAN-CS-72-295  June  1972  28  Pages 


An  Efficient  Parallel  Algorithm  for  the   Solution  of  a  Tridiagonal  Linear 
System  of  Equations 
by  Harold  S.    Stone 

Tridiagonal  linear  systems  of  equations    can  be   solved  on 
conventional   serial  machine   in   a  time  proportional  to  N,   where  N   is 
the  number  of  equations.      The   conventional   algorithms   do  not   lend  them- 
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selves   directly  to  parallel   computation  on   computers   of  the   ILLIAC   IV 
class,    in  the   sense  that  they  appear  to  be   inherently  serial.      An 
efficient  parallel   algorithm  is  presented  in  which  computation  time 
grows   as  log  N.      The  algorithm  is  based  on  recursive   doubling  solutions 
of  linear  recurrence   relations,    and  can  be  used  to   solve  recurrence 
relations  of  all  orders. 

STAN-CS-71-251  December  1971  21   Pages 


NONLINEAR   EQUATIONS 


Error  Bounds,   Based  Upon  Gerschgorin's   Theorems    for  the   Zeros  of  a 

Polynomial 

by  Brian  Thomas   Smith 

Suppose  L  distinct   points   z    ,z0,...z      are   given  in  the   complex 
plane  and  associated  with   each  point   z      is    a  positive   integer  M  ,   the 
'multiplicity'    of  z,  ;    These  points   z    ,  with  their  multiplicities  M  , 
are   supposed  to   approximate  the  zeros   of  a  monic  polynomial   P(z). 
The  object   of  this  thesis    is   to  measure  how  well   a  set   of  such  points 
approximates  the   zeros   of  P(z);   that    is,  we  want  to   determine  how 
close   each  point   z      is  to   its    subset   of  M     zeros. 

TOR  Tech.    Rpt .    9  July    1969  137  Pages  [Ph.D.    Thesis] 


An  Algorithm  for  the  Solution  of  a  Quadratic   Equation  using  Continued 

Fractions 

by  Kishor  Shridharbhai   Trivedi 

This  is   an  effort   to   investigate   representations  of  numbers 
other  than  positional  notation   for   computer  arithmetic.      Using   continued 
fraction  representation  of  numbers,    an   algorithm  to   solve   a  limited 
class  of  quadratics   has  been  developed.      This    algorithm  is   suitable   for 
hardware  implementation   and  is   reasonably  efficient.      Feasibility  of 
constructing  an   arithmetic  unit  with   continued  fraction   representation 
depends   on   discovery  of  many  more   such  useful   algorithms  which   can  share 
the  same  hardware. 

UlUC-R-72-525  June  1972  62   Pages  Program        [M.S.    Thesis] 


A  Method  for  Solving  Polynomial   Equations  by   Continued  Fractions 
by  Amnon  Bracha 

A  method  for  the  approximation  of  all  the  real   roots   of  an 
n-order  polynomial   equation   is    developed.      It    is    assumed  that    intervals 
containing  the  solutions   are  known.      Bilinear  transformations   are  used 
to  approximate  the  solution.      Convergence   is   achieved. 

UIUC-R-72-521  July  1972         20  Pages 
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Kublanovskayas   algorithm  for  solving  the  nonlinear  eigenvalue  problem 
by  Per  Lindstrom 

In  the  present    contribution  an  iterative  algorithm  proposed 
by  B.    N.    Kublanovskaya    [Dokl  Akad  Nauk  188,   lOOU-1005    0-969)    compare 
also  SIAM  Journal  of  Num.   Analysis   7,    532-537   0-970))    for  the  solution 
of  the  nonlinear  eigenvalue  problem  T(A]x  =  0   Is  tested  numerically. 
T(A)   is   an  n*n   real   or  complex  matrix  that    depends   analytically  on 
the  eigenvalue  parameter   A.      It    is   our  aim  to   study  different   possible 
algorithms    for  this  problem  in  order  to   assess  their  practical  value, 
and  this    is   a  part   of  this   investigation. 

UMINF-5.71  August   1971  lh  Pages 


Algorithms    for  the  nonlinear  eigenvalue  problem 
by  Axel  Ruhe 

The   following  nonlinear  eigenvalue  problem  is   studied:      Let 
T(A)  be   an  nxn  matrix,   whose   elements    are   analytical    functions  of  the 
complex  number   A.      Seek  A    and  vectors   x  and  y,    such  that   T(A)x  =  0   and 
y^T(A)   =   0.      Several  algorithms    for  the  numerical  solution  of  this 
problem  are  studied.      These  algorithms    are   extensions  of  algorithms 
for  the  linear   eigenvalue  problem  such   as   inverse   iteration  and  the  QR 
algorithm,    and  algorithms  that   reduce  the  nonlinear  problem  into   a 
sequence  of  linear  problems.      Numerical  tests,  performed  in  order  to 
compare  the   different   algorithms,    are  reported,    and  a  few  numerical 
examples   illustrating  their  behaviour  are   given.      An  ALGOL  program  for 
one  of  the  algorithms   is   given  as    an  appendix. 

UMINF-7-71  August  1971  30  Pages  Program 
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Multi-Dimensional   Quadrature  Formulae 
by  P.    Keast 

Two   classes   of  formulae,    for  numerical   quadrature   in   several 
variables,    are   considered.      Comparisons   are  made  on  the  basis  of 
accuracy  and  efficiency  of  the  various    formulae  used,    and  on  the 
basis   of  ease  of  obtaining  similar   formulae   in  high   dimension.      A 
new  technique   for  obtaining  a  particular  type  of  quadrature   rule   is 
developed.      Several  numerical   examples   are  given. 

TOR  Tech.    Rpt.    1+0  February  1972  78  Pages  Program 
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ORDINARY   DIFFERENTIAL   EQUATIONS 

Studies   in  the  Numerical   Solution  of  Stiff  Ordinary  Differential  Equations 
by  Wayne   En right 

A  class  of  second  derivative  multi-step   formulas   is    developed. 
The   stability  of  these   formulas   is    investigated  and  they  are   shown  to  be 
suitable   for  stiff  equations.      An   implementation  of  a  variable-order, 
variable-step  method  based  on  this   class   of  formulas   is   described  and 
some  numerical   results   are  presented. 

The  reliability  of  methods   suitable   for  stiff  equations    is    also 
investigated.      The   condition  of  a-acceptability  is   introduced  and  it    is 
shown  that    if  a  method  satisfies   this    condition,   then  for  any  member  of  a 
particular  class  of  stiff  systems,   the  numerical    solution  can  be  guaran- 
teed to  lie  on  the  exact   solution  of  a  slightly  perturbed  system. 

Methods   based  on  the   second  derivative   formulas  mentioned  above 
are   shown  to  be  a-acceptable.      Numerical   comparisons  of  methods   in  terms 
of  both  reliability   and  efficiency   are   considered  and  the  background  is 
established  for  a   comprehensive  comparison. 

TOR  Tech.    Rpt .    h6  October  1972  82  Pages  Program  [Ph.D.    Thesis] 


On  the   Identification  of  Multi-Output   Linear  Time-Invariant   and  Periodic 

Dynamic  Systems 

by  Ahmed  H.    Sameh   and  Walter  L.    Heimerdinger 

In  this  paper  we   describe   an  efficient    computational   algorithm 
for  estimating  the  coefficients  of  the   characteristic  polynomial   of  a 
linear  time-invariant  multi-output    dynamic   system,    using  only  output 
observations,    for  qualitative   analysis   of  the  transition  matrix  or  for 
evaluating  its   eigenvalues.      We   also   give   some   computational   results   of 
the   identification  of  those   systems   using  the  Ho-Kalman   approach. 
Furthermore,    an   identification   scheme   for  high-frequency  periodic 
systems  of  unknown  periods   is    described  in   detail. 
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Tchebycheffian  Multistep  Methods    for  Ordinary  Differential   Equations 
by  Tom  Lyche 

The  purpose  of  this  paper  is  to  extend  some  of  the  theory   for 
linear  multistep  methods   to  include   steps ize-dependent    coefficients.      In 
particular  we  treat  the   case  where  we   demand  exact    integration  of  a  given 
set   of  linearly  independent    functions. 
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Numerical  Methods   for  the   Identification  of  Differential  Equations 
by  Raymond  Jonathan  Lermit 

This   dissertation   considers   computational   methods   designed 
to  aid  in  mathematical  model  building.      Specifically,    it   discusses 
methods   of  determining  ordinary  differential   equations   given  their 
solution   in  the   form  of  observed  data. 

Since  the  problem  cannot  be   solved  in  this   generality,    it 
is  necessary  to   supply  equations   containing  arbitrary   functions. 
The  problem  is   then  to   find  these   functions   given  the  solution  of 
the  equations.      In  order  to  be   amenable  to   computer  solution,    a 
discretization   of  the   functions    as    a  linear  sum  of  a  given  ortho- 
normal   set      is   necessary.      The  problem  thus   reduces  to  one   of  find- 
ing a  finite  number  of  parameters. 

The   solution  technique   is  to   find  that   function  which 
produces   a  solution  most    closely   approximating  the  observed  data. 
It    is   thus   a  problem  in  the  minimization  of  nonlinear  functionals 
and  may  be   solved  by   iterative  methods.      Modifications   required  to 
ensure  that  the   functional   is   convex,   thus      guaranteeing  a  global 
minimum  solution,    are   discussed.      Different   algorithms    considered 
for  carrying  out  the  minimization   include  the   generalization  of 
Newton's  method  and  variants   of  it  which   are  more   economical    in 
computer  time,    especially  the  Conjugate   Gradient   method.      All  of 
these  methods   require   derivatives  which   are   derived  automatically 
from  the  original   equation  using  formal   algebraic  manipulation.      The 
different  methods   are   compared  for  rates  of  convergence  and  amount 
of  calculation  required  at   each   iteration. 

Two   examples   are   included.      The  effect   of  introducing 
random  errors   into  the   data  to  simulate  observational   errors,    and 
how  this  may   alter  the   convergence   rates,    is   also   discussed. 
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PARTIAL  DIFFERENTIAL  EQUATIONS 


Finite-Element   Galerkin  Method  for  Mixed  Initial -Boundary  Value 
Problems   in  Elasticity  Theory 
by  Hiroshi   Fujii 

The  Finite-Element   Galerkin  method   (FEG  method)  has   been 
recognized  as   a  very  powerful  tool   for  numerical   solution  of  boundary 
value  problems   of  partial   differential   equations.      In  this   paper,   we 
treat  the   application  of  the  FEG  method  to   vibration  problems    (i.e., 
mixed  initial -boundary  value  problems)    in  linear  elasticity  theory. 
We   introduce   some   special  versions   of  FEG  schemes    for  approximating 
the   solution  of  our  problem,    including  schemes  of  explicit   type   as 
well   as    an   implicit    scheme   and  a   continuous -time   scheme. 
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Automatic   Solutions  of  Partial   Differential   Equations 
by  Leonard  Andrew  Lars en 

The  problem  of  the  present   paper  is  to  look   for  an   automatic 
method  that    can  be   applied  to  some   subset   of  partial   differential   equa- 
tions.     The  primary   value  of  a  program  involving  an  automatic  method 
would  be  to  provide  the  person  who  needs   only  a  few  runs  with   a  parti- 
cular type  of  equation  the  opportunity  to   find  a  solution,   within   accepta- 
ble tolerances,  without   having  to  write   and  debug  a  specialized  program. 
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An  Application  of  Semi- Iterative  and  Second-Degree  Symmetric  Successive 
Overrelaxat  ion   Iterative  Methods 
by   Tran  Phien 

The  goal  of  this   paper  is   to   study  the  effectiveness   of  semi- 
iterative  methods    and  second-degree  methods  based  on  the  symmetric 
successive  overrelaxat  ion  method  for   solving  a  certain  class   of  linear 
systems.      The  linear  systems    correspond  to  the   solution  by   a  finite 
difference  method  of  a  boundary  value  problem  associated  with   a  particular 
self-adjoint  elliptic  partial   differential    equation.      We   shall  use   a  semi- 
iterative  method  based  on  the  symmetric   successive   overrelaxat ion    (SSOR) 
method.      The   resulting  method  is   called  the  SSOR-SI  method.      We  shall 
also   use  the  SSOR-SD  method,    a  stationary   second-degree  method  based  on 
the  SSOR  method. 
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Numerical   Implementation  of  the  Schwarz   Alternating  Procedure   for 
Elliptic  Partial  Differential   Equations 
"by  David  Ross   Stoutemyer 

This  thesis   describes  numerical    implementation  of  the   Schwarz 
and  Neumann   alternating  procedures   for  the  solution  of  the  Laplace-Dirichlet 
problem  on  the  union  of  two   disks,   the   intersection  of  two   disks,    an 
arbitrary  quadrilateral,    and  the  union  of  two   spheres. 

All  of  the   above  examples   lead  to   a  pair  of  coupled  Fredholm 
integral  equations   of  the   second  kind  with   singular  kernels    and  singular 
low-order  derivatives   of  the  solution  at    corners  of  the   region.      These 
singularities   are  overcome  by   a  change  of  variable  together  with  special 
spacing  of  the  abscissas   and  extrapolation  to  the  limit.      These  methods 
are   suitable   for  more  general  boundary   conditions,   more  general  partial 
differential   equations,    and  more  general   geometrical   configurations. 
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Use  of  Fast   Direct   Methods   for  the   Efficient   Numerical   Solution  of 
Non separable  Elliptic  Equations 
by  Paul   Concus   and  Gene  H.    Golub 

We   study   an   iterative  technique   for  the  numerical   solution  of 
strongly   elliptic   equations  of  divergence   form  in  two   dimensions  with 
Dirichlet  boundary   conditions   on   a  rectangle.      The  technique   is  based  on 
the  repeated  solution  by   a  fast   direct   method  of  a  discrete  Helmholtz 
equation  on    a  uniform  rectangular  mesh.      The   problem  is   suitably  scaled 
before   iteration,    and  Chebyshev  acceleration   is    applied  to   improve  conver- 
gence.     We   show  that    convergence   can  be  exceedingly  rapid  and  independent 
of  mesh  size   for   smooth   coefficients.      Extensions   to  other  boundary  condi- 
tions,  other  equations,    and  irregular  mesh  spacings    are   discussed,    and 
the  performance  of  the  technique   is   illustrated  with  numerical   examples. 

STAN-CS- 72-278  April   1972  39  Pages 


Numerical  Solution  of  First-Order  Hyperbolic  Systems   of  Partial 
Differential  Equations 
by  Sunil  K.    Pal 

This  work  develops  two   new   finite-difference  schemes   -  an 
explicit   scheme  and  an   implicit    scheme  -   for  numerical   solution  of  first- 
order  hyperbolic   systems   of  partial   differential   equations   in   any  number 
of  space  variables. 
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INTEGRAL  EQUATIONS 


Methods  for  the  Numerical  Solution  of  Integral  Equations  of  the  Second 

Kind 

by  David  Blair  Coldrick 

A  detailed  analysis   of  the  quadrature  method  is   given   from 
several   points  of  view,    viz.    those  of  Kantorovich,    Buckner,    and  the 
Mysovskih-Brakhage-Anselone   development.      It    is   shown  that  the  technique 
of  deferred  approach  to  the  limit    is   valid  under  fairly  general   circum- 
stances.     An   error  analysis  of  projection    (and  "shifted"  project  ion  )methods 
is   given.      A   class   of  degenerate  kernel  methods    is  proposed,    and 
is   compared  to   closely  related  projection  methods   in  terms   of  error  bounds 
and  ease  of  implementation.      Finally,    a  discussion  of  the   problem  of 
(weak)    singularities   is   presented.      The  application  of  projection  methods 
to   such   equations   is  mentioned.      The  technique  of  product    integration  is 
applied  to  this  problem,    and  a  fairly  general   convergence  statement    is 
established.      For  singularities  the  kernel  k(s,   t)    along  the  line  s   =  t, 
a  generalization  of  Kantorovich's   method  with   a  higher  order  of  conver- 
gence  is   proposed,    and  a  theorem  to  this   effect    is  proved. 
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INTERVAL  AND   SIGNIFICANT  DIGIT  ARITHMETIC 

Implementation  of  Basic  Software   for  Significant   Digit  Arithmetic 
by  Steven  See  Sun  Lai 

This   report    is   concerned  with  the   implementation  of  significant 
digit    arithmetic  using  the  unnormalized  arithmetic   of  Metropolis   and 
Ashenhurst.      One  of  the  primary  objectives  was  to   design  the  basic   soft- 
ware modules   for  inclusion   into  the  OL/2   array  language;   however,   these 
modules   are  written   in   assembly  language   and  therefore   are  adaptable  to 
other  software   systems    for  the   IBM   360/370  machines.      A  discussion  of 
significance  arithmetic,    including  the  nontrivial   problem  of  input/output, 
is  presented.      Examples    are  provided  in  Appendix  A. 
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A  Univac  1108  Program  for  Obtaining  Rigorous   Error  Estimates   for  Approxi- 
mate Solutions  of  Systems  of  Equations 
by  Dennis   Kuba  and  L.    B.    Rail 

A  UNIVAC  1108  computer  program  which  obtains   rigorous   interval 
error  bounds    for  approximate   solutions   of  finite   systems  of  nonlinear 
equations    is   described  in  this   report.      Since  the   coefficients  of  the 
original   system  may  take  on   interval   values,   the   error  bounds   obtained 
include  the  contributions   of  truncation  error  for  Newton's  method,    round-off 
error,    and  possible   errors    in  the  coefficients  of  the   given  system  of  equations 
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OPTIMIZATION 


Large-Scale  Linear  Programming  Using  the   Cholesky  Factorization 
"by  M.    A.    Saunders 

A  variation  of  the  revised  simplex  method  is  proposed  for 
solving  a  standard  linear  programming  problem.      The  method  is   derived 
from  an   algorithm  recently  proposed  by  Gill   and  Murray,    and  is  based 
upon  the   orthogonal   factorization  B  =  LQ  or,    equivalently,   upon  the 
Cholesky  factorization  BB     =  LL     where  B  is  the  usual   square  basis, 
L  is   the  lower  triangular  and  Q  is  orthogonal. 

We  wish  to  retain  the   favorable  numerical  properties  of 
the   orthogonal    factorization,   while   extending  the  work  of  Gill   and 
Murray  to  the   case  of  linear  programs  which  are  both  large  and  sparse. 
The  principal  property  exploited  is  that  the  Cholesky  factor  L   depends 
only  on  which  variables   are   in  the  basis,    and  not   upon  the  order  in 
which  they  happen  to   enter.      A  preliminary  ordering  of  the  rows   of  the 
full   data  matrix  therefore  promises   to  ensure  that   L  will  remain 
sparse  throughout  the  iterations   of  the   simplex  method. 

An   initial    (in-core)   version  of  the   algorithm  has  been   imple- 
mented in  Algol  ¥  on  the  IBM  360/91   and  tested  on  several  medium-scale 
problems    from  industry    (up  to  930   constraints).      While  performance  has 
not  been   especially  good  on  problems  of  high   density,   the  method  does 
appear  to  be  efficient   on  problems  which   are  very  sparse,    and  on 
structured  problems  which  have   either  generalized  upper  bounding, 
block-angular,    or  staircase   form. 

STAN-CS-72-252  January  1972  60  Pages 


Product   Form  of  the  Cholesky  Factorization  for  Large-Scale  Linear 

Programming 

by  Michael  A.    Saunders 

A  variation  of  Gill   and  Murray's   version  of  the  revised  sim- 
plex algorithm  is   proposed,    using  the  Cholesky  factorization  BB1  =   LDL 
where  B  is   the  usual  basis,    D   is   diagonal   and  L   is  unit   lower  triangular. 
It    is   shown  that    during  change  of  basis   L  may  be   updated  in  product   form. 
As  with   standard  methods   using  the  product   form  of  inverse,   this   allows 
use   of  sequential   storage   devices   for  accumulating  updates   to  L.      In 
addition  the   favorable  numerical  properties  of  Gill   and  Murray's   algorithm 
are  retained. 

Close   attention   is   given  to   efficient   out-of-core   implementation. 
In  the   case  of  large-scale  block-angular  problems,   the  updates  to  L  will 
remain  very  sparse   for  all   iterations. 
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The  Differentiation  of  Pseudoin verses   and  Nonlinear  Least  Squares 
Problems  Whose  Variables  Separate 
by  G.    H.    Golub   and  V.    Pereyra 

For  given   data   (t . ,  y.),    i  =  1,    .  . . ,  m,  we  consider  the 
least   squares   fit-  of  nonlinear  models  of  the  form 

n 
F(a,    a;   t)   =      I     g.(a)    0.(a;   t),    a   eflf,    a   eR,  . 
~  j=l     J    ~        J 

For  this   purpose  we  study  the  minimization  of  the  nonlinear  functional 

m  2 

r(a,    a.)   =     I      (y.    -  F(a,    a,   t.))    . 
i=l       1  1 

It    is   shown  that   by  defining  the  matrix{$(a) }.     .    =    <j>.(a;  t.)    and  the 

+      ~      i,J         \]    ~        i 

modified  functional  r    (a)   =  ||y  -   $(a)$      (a)   y||^    ,    it    is  possible  to 

optimize   first  with  respect   to  the  Parameters   a,    and  then  to  obtain, a 
posteriori,   the  optimal   parameters   a.      The  matrix  $+(g)    is  the  Moore- 
Penrose   generalized  inverse  of  $(ot),    and  we   develop   formulas   for  its 
Frechet    derivative  under  the  hypothesis  that   $(a)    is   of  constant    (though 
not   necessarily  full)   rank.      From  these   formulas  we  readily  obtain  the 
derivatives   of  the   orthogonal   projectors   associated  with   $(a),    and  also 
that   of  the   functional   r9(a).      Detailed  algorithms   are  presented  which 
make  extensive  use  of  well-known   reliable  linear  least    squares  techniques, 
and  numerical   results   and  comparisons    are  given. 
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The  method  of  conjugate   gradients   used  in  iverse   iteration 
by  Axel  Ruhe  and  Torbjorn  Wiberg 

An   algorithm  is   devised  that    improves   an  eigenvector  approxi- 
mation  corresponding  to  the  largest    (or  smallest)    eigenvalue  of  a  large 
and  sparse  svmmetric  matrix.      It   solves  the  linear  systems  that   arise 
in   inverse      iteration  by  means  of  the   c-g   algorithm.      Stopping  criteria 
are   developed  which  insure   an   accurate  result,    and  in  many   cases    give 
convergence   after  a  small  number  of  c-g  steps. 

UMINF-22.72  April   1972  21  Pages 


Partial  Analysis   of  an  Algorithm  for  Minimizing  Functions  Based  on  a 
Homogenous  Model 
by  James  ¥.    Daniel 

Jacobson   and  Oksman  recently  presented  a  method  which  they 
claimed  minimized  a  homogeneous    function  of  arbitrary   degree   in   finitely 
many  steps    and  which  was   globally   convergent    for  a  large   class   of  general 
functions.      We  point   out  here  errors   in  the  proofs   of  these  results, 
suggest    some  modifications    in  the  method  so  that  the  global    convergence 
result    is   valid,    and  we  give   some   examples   to   show  the   computational 
effects   of  these   changes. 
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Convergent   Step  Sizes   for  Curvilinear-Path  Methods  of  Minimization 
by  James  W.    Daniel 

When  one  seeks   a  point  minimizing  a  function  f  over  a  set 

C — which  may  be  the  whole   space — one  often  moves    from  one  approximate 

solution  x     to   another  x     .,    =  x     +tp     by  searching  along  a  ray 
n  n+1  n  n  n 

x     +  tp   ,   where  the  step-size  t     must  be  chosen.      More  generally  one 
n  n  n 

can   consider  moving  along  a  curvilinear  path  x     +  p    (x   ,t ) .      The 

n         n     n 

present   paper  shows  how  the  usual   step-size  algorithms   can  be  used 
with  this   general   approach. 
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Nonlinear  Least   Squares   Without   Derivatives:      An  Application  of  the 
QR  Matrix  Decomposition 
by  Richard  H.    Bartels 

A  new  algorithm  is  proposed  for  minimizing  sums   of  squares 
of  nonlinear  functions   of  several  variables  without  the  use  of  deriva- 
tives.     The   algorithm  is   constructed  using  a  multivariate  secant 
technique    (Broyden  rank-one  method)   to   approximate  the  Jacobian  of 
the   summands,    and  this   Jacobian  is   used  to   drive  the  Levengerg-Marquardt 
iteration.      The  Jacobian  is  kept   and  updated  in   QR  decomposition,    and 
the  linear  least   squares   subproglems  which   arise   at   each  iteration   cycle 
are   solved  by  Golub's  method.      A  FORTRAN  program  and  test     results    are 
offered.      It   is   shown  that  the  algorithm  accomodates   linear  equality  con- 
straints with  ease.      Remarks    are  made   concerning  the  treatment   of  rank 
deficiencies   in  Golub's  method  with   respect  to   algorithms  of  the  type 
Presented. 
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Global   Convergence   for  Newton  Methods   in  Mathematical  Programming 
by  James  W.    Daniel 

In  constrained  optimization  problems    in  mathematical  program- 
ming one  wants   to  minimize   a  functional    f(x)   over  a  given   set   C.      If, 
at    an   approximate  solution   x   ,    one  replaces    f(x)   by  its   Taylor  series 
expansion  through  quadratic  terms   at   x     and  denotes  by  x  the  minimiz- 

ing point    for  this   over  C,    one  has    a  direct   analogue  of  Newton's  method, 
The  local   convergence  of  this  has  been  previously  analyzed;   here  we 
give  global   convergence   results    for  this   and  the  similar  algorithm  in 
which  the   constraint   set   C   is    also  linearized  at   each   step. 
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Newton's  Method  for  Nonlinear   Inequalities 
by  James  W.    Daniel 

A  Newton-type  algorithm  has  "been  presented  elsewhere  for 
solving  non-linear  inequalities  of  the   form  f(x)   ^  0,    and  quadratic 
convergence  has  been  proved  under  very  strong  hypotheses.      In  this 
paper  we   show  that  the   same  results  hold  under  a  considerable  weaken- 
ing of  the  hypotheses. 
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GRAPH  ALGORITHMS 


Graph   Isomorphism 

by  Derek  Gordon   Corneil 

A  procedure   for  determining  whether  two   graphs   are   isomorphic 
is   described.      During  the  procedure,    from  any  given  graph  two   graphs, 
the  representative  graph   and  the   reordered  graph,    are   derived. 

The  time  required  to   determine  both   derived  graphs   depends   on 
a  power  of  n,   the  order  of  the   given  graph.      This   power  is    a  function 
of  an   adjacency  property  known   as  the   strong  regularity  of  the   given 
graph.      For  graphs  that    do   not   contain   a  strongly  regular  transitive 
subgraph,   the  power  is,    at   worst,    five. 

All   given  timing  estimates   for  graphs  that    do  not    contain   a 
strongly  regular  transitive   subgraph   are   confirmed.      The   algorithm  has 
been   programmed  and  in  the   implemented  version  one   of  the   following  four 
messages  will   come  out:       (l)    The   graphs    are   isomorphic.      (2)    The   graphs 
are  not    isomorphic.       (3)   The  representative   graphs   are   identical;   the 
reordered  graphs    are  not    identical;   hence,   these   graphs   form  a  counter- 
example to  the   conjecture.      (k)   The   graphs   contain   a  3-strongly  regular 
subgraph. 
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An  Algorithm  to  Determine  the   Chromatic  Number  of  a  Graph 
by  Barry  Graham 

A  heuristic  algorithm  for  the   determination  of  the   chromatic 
number  of  a  finite   graph  is   presented.      This    algorithm  is  based  on   Zykov's 
theorem  for   chromatic  polynomials   and  extensive  empirical  tests   show  that 
it    is   the  best    algorithm  available.      Christofides '    algorithm  for  the 
determination  of  chromatic  number  is    described  and  is   used  in  the   compari- 
son tests.      An  Algol- W  coding  of  both   algorithms    is   included  in  the   appendix. 
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Algorithms   for  Finding  Cliques  of  a  Graph 
by  Gordon  D.    Mulligan 

Various  methods   for  determining  cliques   in  undirected  graphs 
are  presented  and  analyzed.      Testing  schemes  to   compare  the  methods  on 
graphs  with  a  maximum  number  of  cliques   and  on  graphs  that    attempt   to 
represent   some   applications  are  described. 

A  theorem  states  that  the  maximum  number  of  cliques   in   a  graph 
increases   exponentially  as  the  number  of  vertices   increases.      For  this 
reason  heuristics    are  employed  in   clique   finding  algorithms  to   decrease 
the  amount   of  search   and  an  attempt    is  made  to  measure  the   algorithms   as 
a  function  of  the  number  of  cliques. 

Three   algorithms   are   formally  defined  and  tested:      the  Bier- 
stone   algorithm,   the  Bron-Kerbosch   algorithm,    and  the  Corneil   algorithm. 
The  Bron-Kerbosch   algorithm  is   judged  to  be  the  best. 
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Spectra  of  Finite  Graphs 

by   L.    Coll at z,    Hamburg  and  U.    Sinogowitz,    Darmstadt 

Various   problems   encountered  in  practical   applications,    such 
as   approximations   of  characteristic   frequencies  of  a  membrane  with 
fixed  perimeter  and  given  area,    or  calculations  of  air  vibrations    in   a 
bounded  volume  using  difference  methods,    gave   rise  to  the   consideration 
of  graphs   in   general. 

In  this   paper   finite    (in  a  subsequent   paper  certain  types   of 

infinite)    graphs   are  examined  and  in  particular   some   conclusions   drawn 

with  regard  to  their  relationship  to  the  theory  of  indecomposable,   non- 
negative  matrices. 
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Graphs,    Groups    and  Matrices 
by  Abbe  Mowshowitz 

Our  object  here  is  to  exploit  the  connection  between  the 
adjacency  matrix  of  a  graph  and  its  automorphism  group  in  order  to 
determine  the  latter. 

In  what   follows,   we  will   take  advantage  of  the   fact   that  the 
adjacency  matrix  of  a  graph  is    a   (0,    l)-matrix   and  thus    can  be  regarded 
as   a  matrix  over  GF(2). 

UBC   Tech.    Rpt.    71-OU  October  1971  1^   Pages 
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Combinatorial  Solutions   to  Partitioning  Problems 
by  J.    A.    Lukes 

In  this   dissertation  we   describe   algorithms  that   use   graph 
properties   and  dynamic  programming  techniques  to  generate  the  optimal 
partition  of  an   arbitrary  graph.      In  particular,   let    G  be   a  graph 
with  weighted  nodes   and  weighted  edges.      We   consider  algorithms  that 
solve  the  problem  of  partitioning  G  into   clusters  of  nodes   such  that 
the  sum  of  the  node  weights   in  any  cluster  does  not   exceed  a  given  maxi- 
mum W  and  the  weights   of  the   intercluster  edges   are  minimized.      An 
interesting  application  of  such   an  algorithm  is  the   assignment   of  a 
program's    subroutines    and  data  to  pages   in   a  paged  memory  system  so 
as   to  minimize  paging  faults. 

A  very  efficient   variation  of  the  general   algorithm  results 
if  the   graph  to  be  partitioned  is   a  tree.      We   show  that   trees    can  be 
partitioned  in  a  time  proportional  to  the  number  of  nodes   in  the  graph. 

STAN-CS-72-293  May  1972  120  Pages  [Ph.D.    Thesis] 


Chromatic  Automorphisms  of  Graphs 
by  V.    Chvatal   and  J.    Sichler 

The   coloring  group  and  the   full    automorphism  group  of  an 
n-chromatic  graph   are   independent    if  and  only  if  n   is    an   integer  >    3< 

STAN-CS-72-273  March  1972  11  Pages 
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